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I.  INTRODUCTION 


Accurately  characterizing  particulate  mixtures  is  critical  to  a  variety  of  scientific  fields 
including  heterogeneous  materials  like  solid  rocket  propellants  [1,  2],  granular  materials  like 
sand  [3],  low-temperature  phases  of  matter  [4],  and  biological  applications  such  as  protein 
configurations  [5].  Processing  plays  a  critical  role  in  material  science  also.  It  is  known 
that  material  manufacturing  steps  lead  to  anisotropic  morphologies.  For  example,  in  solid 
propellants  an  uneven  burning  profile  is  observed  due  to  particle  alignment  across  the  propel¬ 
lant  grain,  which  is  attributed  to  the  casting  procedure  [6].  Moreover,  determining  optimal 
packings  of  various  shapes  is  a  fundamental  mathematical  problem  that  has  been  studied  for 
centuries  [7,  8].  Often,  packing  algorithms  are  used  to  study  how  particles  arrange  and  are 
the  basis  for  analyzing  the  effect  of  the  morphology  on  physical  phenomena  [9-11].  However, 
Caulkin  et  al.  [12]  showed  that  one  must  consider  physical  interaction  forces  due  to  material 
formulation  steps  in  these  algorithms  to  model  real  microstructures  accurately.  This  is  a 
difficult  task  computationally.  Therefore,  micro-computer  tomography  (micro-CT)  has  be¬ 
come  a  popular  method  for  obtaining  a  description  of  real  microstructures  including  solid 
propellants  [13,  14],  glass  beads  [15,  16],  and  Fontainebleau  sandstone  [17],  just  to  name 
a  few.  Moreover,  tomography  characterization  methods  can  be  used  to  validate  packing 
algorithms. 

For  analyzing  various  material  systems,  statistical  characterization  methods  have  been 
used  to  describe  these  complicated  microstructures.  Bernal  [18]  motivated  the  importance 
of  higher-order  statistics  using  the  radial  distribution  function  to  describe  the  structures 
of  liquids.  Frisch  and  Stillinger  [19]  introduced  n-point  probability  functions  to  describe 
radiation  scattering  in  packs  of  monodisperse  spheres.  Note  that  n-point  probability  func¬ 
tions  can  be  used  to  statistically  characterize  a  general  class  of  microstructures  including 
porous  materials  and  particulate  packs  of  arbitrarily  shaped  inclusions.  The  significance  of 
the  statistical  description  has  also  been  shown  in  several  other  fields  of  physics  including 
condensed  matter  physics  with  applications  to  disordered  materials  through  non-Gaussian 
noise  [20]  and  to  cosmic  radiation  using  Minkowski  functionals  [21],  While  these  statis¬ 
tical  descriptors  have  been  used  for  decades,  accurately  obtaining  higher  order  statistical 
information  of  real  systems  in  three  dimensions  has  proved  difficult  to  this  day.  Coker  and 
Torquato  [22]  and  Yeong  and  Torquato  [23]  computed  two-point,  probabilities  functions  of  a 
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two-phase  material,  yet  considered  only  two-dimensional  slices  and  assumed  rotational  and 
translational  invariance  of  the  functions.  Fullwood  et  al.  [24]  computed  n-point  probability 
functions  of  three-dimensional  two-phase  materials,  but  their  computations  were  limited  to 
relatively  small  systems.  In  [15,  25],  radial  distribution  functions  are  used  to  statistically 
characterize  the  microstructures,  which  are  limited  to  packs  of  spheres.  As  multi-phase  poly- 
disperse  systems  are  considered  in  three  dimensions,  large  packs  must  be  used,  resulting  in 
even  larger  datasets.  In  our  prior  work  [16],  one-,  two-,  and  three-point  probability  functions 
were  calculated  for  a  three-dimensional  system  of  polydisperse  spheres  from  micro-CT  data. 
Mecke  and  Arns  [26]  alternatively  used  Minkowski  functionals  as  statistical  descriptors  to 
describe  porous  materials,  such  as  soils  and  sedimentary  rocks.  The  Minkowski  function¬ 
als  are  better  suited  for  two-phase  materials  characterized  by  a  Boolean  model  in  which  a 
material  skeleton  is  represented  with  overlapping  inclusions. 

Linking  the  morphology  to  the  overall  material  response  is  a  longstanding  problem.  The 
seminal  multiscale  techniques  proposed  by  Hashin  and  Shtrikman  [27],  Willis  [28],  Bensous- 
san  et  al.  [29],  Hill  [30],  and  Torquato  [31]  deserve  attention,  and  many  others  have  also 
contributed  [32-35].  As  computationally  and  tomographically  derived  packs  are  often  too 
large  for  practical  numerical  simulations  of  nonlinear  processes  [1,  36],  work  has  been  ded¬ 
icated  to  reconstructing  statistically  equivalent  unit  cells  [16,  37].  In  contrast  to  numerical 
simulations  [1,  36],  analytical  and  semi-analytical  techniques  are  an  alternative  means  of 
computing  the  overall  response  of  heterogeneous  materials.  The  work  of  Willis  is  of  par¬ 
ticular  importance  as  it  established  a  direct  link  between  the  statistical  description  and 
properties  [28].  However,  in  [28]  results  were  only  obtained  for  simple  morphologies  with 
statistically  isotropic  or  ellipsoidal  symmetric  statistics.  This  limitation  is  due  to  challenges 
in  integrating  complex  integral  kernels  that  are  products  of  the  second  derivative  of  Green’s 
function  and  two-point  probability  functions.  Drugan  and  Willis  [38],  when  formulating  a 
non-local  model,  were  limited  to  obtaining  overall  properties  for  two-phase  systems  with 
an  isotropic  distribution  of  monomodal  spheres.  In  the  work  of  Torquato  and  Sen  [39], 
only  statistically  isotropic  microstructures  and  those  composed  of  aligned  arbitrarily  shaped 
inclusions  were  discussed.  Ponte  Castaneda  [40]  among  others  [38,  41]  have  proposed  a 
third  order  model  for  nonlinear  materials,  yet  once  again  only  simple  microstructures  were 
considered  due  to  the  difficulty  in  obtaining  real  probability  functions  and  integrating  over 
them. 
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In  this  work,  voxel  datsets  consisting  of  polydisperse  particles  of  multiple  shapes  are 
studied.  Using  our  SHAPE3D  software  package,  particles  in  a  voxel  dataset  are  modeled  as 
ellipsoids  by  minimizing  differences  in  geometric  quantities  between  the  voxel  and  Euclidean, 
M3,  representations.  This  description  in  continuous  Euclidean  space  can  simplify  analysis 
of  such  systems,  improve  their  understanding,  and  aid  in  the  development  of  new  material 
formulations.  Moreover,  an  idealized  mapping  can  reduce  the  size  of  the  dataset.  This  re¬ 
duction  becomes  increasingly  important  as  more  complicated  systems  are  studied,  and  larger 
voxel  datasets  must  be  considered  in  order  to  resolve  all  important  morphological  features. 
Note  that  state  of  the  art  in  micro-CT  analysis  of  highly  packed  particulate  systems  has  been 
limited,  to  the  best  of  our  knowledge,  to  nearly  monodisperse  packs  [12,  15,  25],  and  these 
works  do  not  discuss  in  detail  image  processing  steps.  However,  image  processing  steps  can 
pollute  mechanical  and  transport  properties  predictions  that  are  obtained  from  improperly 
segmented/processed  images.  When  considering  polydisperse  packs  of  various  inclusions,  a 
traditional  image  processing  pipeline  does  not  properly  identify  individual  particles.  Thus, 
a  novel  image  processing  pipeline  is  presented. 

Next,  we  compute  overall  anisotropic  properties  without  traditional  assumptions  on  the 
complexity  of  the  microstructure  (statistically  isotropic,  aligned  inclusions).  We  show  that 
for  systems  of  polydisperse  unaligned  inclusions,  the  importance  of  eliminating  these  as¬ 
sumptions  is  critical.  The  n-point  probability  functions  are  computed  directly  from  three- 
dimensional  micro-CT  datasets  using  our  code,  STAT3D  [16].  We  consider  polydisperse 
systems  and  present  that  statistical  functions  are  not  isotropic  nor  possess  ellipsoidal  or 
any  other  material  symmetry.  Thus,  the  closed-form  solutions  of  overall  properties  are  un¬ 
achievable.  Therefore,  based  on  our  prior  work  [42],  the  complex  integrals  containing  the 
n-point  probability  functions  and  the  second  derivative  of  the  Green’s  function  are  com¬ 
puted  numerically.  Various  compositions  of  polydisperse  mixtures  of  spheres  and  ellipsoids 
are  considered,  and  overall  anisotropic  thermal  conductivity  tensors  are  predicted. 


II.  SAMPLE  PREPARATION,  MICRO-CT  SCANNING,  AND  PARTICLE  CHAR¬ 
ACTERIZATION 

This  section  describes  the  methods  to  mix  and  pack  polydisperse  systems,  to  acquire 
data  using  micro-CT,  to  process  datasets  using  image  processing  algorithms,  and  to  model 
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particles  as  ideal  shapes.  In  this  work,  a  surrogate  system  of  rice  grains  and  mustard  seeds 
is  used  for  their  ellipsoidal  and  spherical  shape,  respectively.  In  our  discussion,  we  label  the 
mustard  seeds  as  spheres  due  to  their  low  eccentricities.  However,  all  particles  are  modeled 
as  ellipsoids. 

A.  Packing  of  Polydisperse  Particulate  Mixtures 

Two  packing  methods  were  developed  with  previous  studies  in  mind  [15,  43,  44],  The 
first  method  is  used  for  packing  homogeneous  polydisperse  mixtures  of  ellipsoids  and  spheres 
while  limiting  boundary  effects.  Packs  mixed  using  this  method  will  be  referred  to  as  ran¬ 
domized  packs.  The  second  method  is  used  for  packing  ordered  systems  of  particles,  with 
the  goal  of  obtaining  aligned  inclusions.  Packs  created  with  aligned  particles  are  referred  to 
as  semi-ordered  packs. 


FIG.  1.  The  scanning  container  filled  with  50%  weight  of  spheres  to  ellipsoids.  Note  the  hemi¬ 
spherical  beads  on  the  boundary  of  the  container  used  to  minimize  boundary  effects. 

First,  the  preparation  of  randomized  packs  of  polydisperse  systems  is  described.  The 
container,  as  seen  in  Figure  1,  was  designed  to  eliminate  boundary  effects  and  fill  the  entire 
viewing  area  within  the  micro-CT  scanner.  The  container  has  a  diameter  of  62  mm  and 
height  of  65  mm,  and  hemispherical  beads  (3  mm  radius)  are  attached  to  the  boundary  in 
order  to  minimize  boundary  effects.  The  randomized  packs  are  prepared  in  the  following 
steps:  I)  Amount  of  each  constituent  (e.g.  spheres,  ellipsoids,  etc.)  is  weighed  with  a  high 
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precision  scale  with  0.001  g  accuracy.  II)  Constituents  are  mixed  together  in  a  large  bowl. 
Ill)  Mixture  is  poured  into  scanning  container.  IV)  Scanning  container  is  shaken  in  all 
directions.  V)  Container  is  tapped  allowing  particles  to  settle  and  reach  maximum  volume 
fraction.  Let  us  note  that  steps  II  and  IV  are  repeated  until  the  mixtures  are  visually 
homogeneous. 

To  prepare  semi-ordered  packs  of  aligned  particles,  a  separate  sample  preparation  strategy 
is  used.  The  same  sized  cylindrical  container  was  used,  but  no  hemispherical  beads  were 
attached  to  the  boundary  of  the  container  in  order  to  prevent  the  irregular  boundary  from 
randomizing  the  sample.  The  semi-ordered  packs  are  prepared  in  the  following  steps:  I)  A 
thin  layer  of  particles  is  poured  into  container.  II)  On  a  flat  surface,  container  is  oscillated 
back  and  forth  twenty  times  in  the  same  direction.  Ill)  Tap  container  five  times.  These 
steps  are  applied  to  thin  layers  so  that  the  particles  on  top  of  the  pack  can  move  freely.  We 
repeat  these  steps  until  the  container  is  full.  The  oscillation  in  the  same  direction  allows 
particles  to  align,  and  the  tapping  is  again  done  to  reach  a  maximum  volume  fraction. 


B.  Dataset  Acquisition 

Once  the  packs  are  prepared,  all  samples  are  scanned  using  a  Skyscan  1172  micro-CT 
machine.  This  particular  scanner  has  the  ability  to  produce  datasets  that  capture  features 
that  are  ~  0.7  /an.  A  convergence  study  was  done  on  a  pack  of  polydisperse  spheres  to  show 
that  the  particle  volume  fraction  of  a  sample  is  maintained  for  resolutions  of  69.4,  34.7,  17.4, 
and  8.7  //m  per  pixel.  The  standard  deviation  in  volume  fraction  for  all  of  these  scans  was 
0.0055,  which  is  less  than  1%  of  the  mean  volume  fraction.  A  resolution  of  69.4  /irn  per  pixel 
was  selected  for  all  compositions  described  in  this  work.  All  datasets  are  on  the  order  of 
108  voxels  (voxel  =  3D  pixel).  The  average  diameter  of  the  spherical  particles,  the  smallest 
sized  inclusion  in  our  mixtures,  is  ~  2  mm.  With  a  resolution  of  69.4  /irn,  this  is  ~30  voxels 
in  diameter,  which  is  a  good  resolution  for  the  statistical  characterization  that  follows. 

The  three-dimensional  voxel  dataset  of  the  50%  weight  fraction  of  spherical  to  ellipsoidal 
particles  is  shown  in  Figure  2.  This  dataset  will  be  used  to  describe  steps  for  characterizeing 
the  microstructure.  This  cylindrical  dataset  has  a  diameter  of  D  =  57.67  mm  (831  pixels) 
and  a  height  of  H  =  48.51  mm  (699  pixels),  which  corresponds  to  3.79  x  108  voxels. 
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FIG.  2.  The  surrogate  system  of  ellipsoids  and  spheres  with  50%  weight  fraction  of  constituents. 
Observe  the  hemispherical  beads  at  the  boundary  that  are  used  to  randomize  the  packs  and  limit 
boundary  effects.  The  brightest  particles  correspond  to  the  ellipsoidal  particles,  while  the  darker 
ones  correspond  to  spherical  ones. 

C.  Image  Processing 

Once  the  three-dimensional  voxel  dataset  is  obtained,  image  processing  algorithms  are 
used  to  identify  particles,  and  these  voxel  particles  are  then  modeled  as  idealized  shapes. 
The  purpose  of  characterizing  particles  as  idealized  shapes  is  to  understand  the  influence 
of  different  particle  types  on  the  macroscopic  behavior  and  to  provide  a  description  of  the 
microstructure  in  continuous  Euclidean  space,  M3.  This  compact  representation  can  be  ben¬ 
eficial  for  design  of  new  material  formulations,  for  example.  Image  processing  routines  for 
identifying  particles  in  a  voxel  dataset  are  well  established  [45].  For  a  system  of  monodis- 
perse,  convex-shaped  objects  with  no  hollow  regions,  a  typical  image  processing  pipeline 
is  to  eliminate  noise  in  the  dataset  by  using  an  edge-preserving  smoothing  algorithm,  to 
identify  the  particulate  material  by  thresholding,  and  to  segment  the  dataset  such  that  indi¬ 
vidual  particles  are  represented  by  connected  groups  of  voxels.  Two  common  segmentation 
algorithms  are  watershed-based  segmentation  and  the  opening  algorithm  (erosion  followed 
by  dilation).  Unfortunately,  a  watershed-based  algorithm  is  limited  as  it  only  properly  sepa¬ 
rates  datasets  of  mono-disperse  convex-shaped  objects,  while  the  opening  algorithm  distorts 
shapes  for  non-spherical  objects.  In  general,  all  segmentation  algorithms  lead  to  some  vol¬ 
ume  removal.  For  polydisperse  systems  with  non-convex  features,  these  standard  image 
processing  steps  will  result  in  improperly  segmented  datasets  [45].  As  real  materials  consist 
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of  irregularities,  such  as  voids,  non-convex  surfaces,  varying  densities,  etc.,  this  traditional 
processing  pipeline  must  be  enhanced  to  accurately  capture  real  microstructures. 


(a)  (b)  (c) 

FIG.  3.  Illustration  of  image  processing  steps  i-iii. 


In  order  to  limit  and  control  the  errors  in  the  image  processing  pipeline,  we  propose  a 
novel  strategy  based  on  interplay  between  voxel  analysis  and  analysis  in  M3.  The  following 
image  processing  steps  were  developed  using  Visage  Imaging  Inc.’s  Amira  [46]  in  conjunction 
with  Skyscan’s  CTAn  [47]  and  the  SHAPE3D  software  developed  in  this  work:  i)  The  voxel 
dataset  is  smoothed  using  the  Kuwahara  edge-preserving  median  filter  algorithm  [48]  in 
order  to  eliminate  noise.  A  resulting  two-dimensional  slice  of  the  dataset  after  smoothing  is 
presented  in  Figure  3(a).  ii )  The  particulate  phases  are  thresholded  based  on  the  grayscale 
distribution  of  particles  (Figure  3(b)).  In  this  distribution,  the  three  peaks  on  the  right 
are  related  to  the  average  grayscale  values  of  the  ellipsoids,  spheres,  and  the  container 
and  the  left  peak  is  considered  void  space.  The  threshold  cutoff  (red  dashed  line),  which 
distinguishes  between  void  space  and  material,  is  chosen  to  be  the  local  minimum  value 
between  the  two  leftmost  peaks  and  is  60  for  this  dataset,  in )  Voids  within  the  particles 
are  filled.  The  resulting  dataset  slice  after  this  step  is  presented  in  Figure  3(c).  Notice 
how  many  of  these  particles  are  joined,  iv)  A  watershed  segmentation  algorithm  is  used  to 
separate  individual  particles,  v)  To  eliminate  the  boundary  effects  due  to  the  hemispherical 
beads  that  are  fixed  to  the  scanning  container,  a  volume  of  interest  (VOI)  is  defined.  The 
VOI  (the  circular  blue  region  in  Figure  4(a))  is  chosen  such  that  the  boundary  effects  on  the 
statistical  characteristics  are  minimized/eliminated.  The  VOI  is  determined  by  decreasing 
the  cylindrical  volume’s  diameter  and  height  by  increments  of  10  voxels  until  the  particle 


volume  fraction  saturates  to  within  1%.  The  resulting  VOI  has  a  diameter  of  -Dyoi  = 
46.57  mm  and  a  height  of  i/yoi  =  37.4  mm  (63.70  cm3),  vi)  Geometric  quantities  of 
individual  particles  including  volume,  surface  area,  moments  of  inertia,  the  mean  grayscale 
value  and  centroid  coordinates  are  computed,  vii )  It  is  determined  which  particles  are  still 


(a)  (b) 

FIG.  4.  Illustration  of  image  processing  steps  iv-vii.  (a)  A  dataset  slice  after  image  segmentation  in 
step  iv.  (b)  A  section  of  (a)  that  illustrates  particles  remaining  connected  after  step  vi  (highlighted 
in  orange). 

connected  due  to  the  inefficiency  of  the  segmentation  algorithm,  and  these  particles  are 
marked  (Figure  4(b)).  To  determine  which  particles  were  not  separated  in  step  iv,  one  can 
analyze  certain  geometric  quantities  of  the  voxel  particles  to  see  if  they  match  the  parameters 
of  the  shapes  expected  in  the  sample,  i.e.  if  some  voxel  particles  are  larger  than  the  known 
distribution  of  particles,  they  should  be  segmented  further.  Because  this  system  is  composed 
of  particles  closely  resembling  ellipsoids,  a  voxel  particle  not  matching  an  ellipsoid  is  marked 
for  additional  image  segmentation  steps.  The  method  for  modeling  a  voxel  particle  as  an 
ellipsoid  using  our  SHAPE3D  package  is  described  in  subsection  II D.  viii )  A  combination  of 
the  opening  algorithm  with  watershed  segmentation  is  used  to  separate  the  marked  particles 
identified  in  step  vii  using  SHAPE3D  (see  subsection  II D).  ix )  For  newly  separated  particles, 
geometric  quantities  are  computed,  x)  The  phase  of  each  particle  is  determined  based  on 
the  mean  grayscale  value.  The  average  grayscale  value  of  each  particle  is  described  in 
Figures  5(b)  and  5(c),  where  d,  is  the  representative  sphere  diameter.  Considering  the  local 
minimum  in  Figure  5(b)  and  the  distribution  of  particle  sizes  in  Figure  5(c),  the  average 
grayscale  cutoff  value  for  the  50%  weight  fraction  of  spheres  is  136. 
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FIG.  5.  Illustration  of  image  processing  step  x. 


Figure  6  shows  the  resulting  dataset  after  image  processing,  where  the  spherical  and 
ellipsoidal  particles  are  colored  blue  and  green,  respectively.  Packs  of  voxel  particles  will  be 
referred  to  as  the  voxel  packs.  An  error  measure  is  introduced  to  quantify  the  volume  loss 
due  to  the  image  segmentation  steps, 

o  <  eIP  =  |Cp’  ~BP'  1  X  100  <  100  [%],  (1) 

Cp 

where  cp,B  and  cp,A  are  the  total  particle  volume  fractions  of  the  voxel  pack  before  and 
after  image  segmentation  steps,  respectively.  For  the  50%  weight  spherical  particle  pack, 
cp,B  =  0.6770  and  cp,A  =  0.6594,  which  represents  a  volume  loss  of  Sip  =  2.60%  as  a  result 
of  the  image  segmentation  algorithms  used.  Note  that  the  respective  phases  (mustard  and 
rice)  can  only  be  identified  after  all  image  segmentation  steps,  and  cve  and  cvs  will  refer  to  the 
volume  fractions  of  ellipsoids  and  spheres  in  the  voxel  pack  (cp,A  =  cve  +c"),  respectively.  For 
this  pack,  cve  =  0.2701  and  cvs  =  0.3893.  Relevant  information  including  the  volume  fractions 
and  other  morphological  characteristics  of  all  mixtures  studied  are  presented  at  the  end  of 
this  section  (see  Tables  I  and  II). 


D.  Modeling  Voxel  Pack  as  Pack  of  Idealized  Shapes 

Taking  into  consideration  knowledge  about  the  microstructure,  all  voxel  particles  are 
modeled  as  ellipsoids  in  M3.  The  voxel  particles  are  fitted  to  idealized  shapes  by  optimizing 
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FIG.  6.  Voxel  pack  of  50%  weight  of  constituent  mixture.  The  voxel  pack  is  the  voxel  dataset  after 
image  processing. 


two  fitness  functions  that  best  characterize  a  particle’s  volume,  surface  area,  and  orientation. 
The  first  function  minimized  is 


niW="  IK, I  +(1^>  \S«\  '  (2) 

where  Vct  and  Sct  are  the  volume  and  surface  area  calculated  for  a  particle  in  the  voxel  pack, 
Vpa  and  Spa  are  the  volume  and  surface  area  of  the  ideal  shape  being  optimized,  and  a;  is  a 
weight  factor  between  a  particle’s  volume  to  surface  area.  The  weight  factor  was  chosen  to 
be  oj  =  0.8  in  our  work.  The  scaling  factor  k  determines  the  scale  of  the  final  semiaxes  as 


a 

a* 

b 

=  k 

b* 

L  C  , 

L  C*  , 

(3) 


where  a,  b,  and  c  are  the  final  semiaxes  of  the  ideal  shape  with  a  >  b  >  c.  The  semiaxes  a*, 
b*,  and  c*  are  computed  from  the  inertia  tensor  of  the  voxel  particle,  Ict,  assuming  that  the 
particle  is  an  ellipsoid.  The  use  of  a  scaling  factor  k  ensures  that  the  ratios  between  a*,  b*, 
and  c*  and  a,  b,  and  c  are  preserved  so  that  the  idealized  particle  better  mimics  the  overall 
shape. 

A  second  fitness  function,  n2(a,/3,7),  is  defined  to  adjust  the  orientation  of  the  idealized 
shape  to  best  mimic  the  voxel  particle  as 


n2  (a,  (3, 7) 


\R(a,/3,-f)  IpaR  (a, /3,^) 


-1 


lct\ 


Lct\ If 


(4) 
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where  Ipa  is  the  inertia  tensor  of  the  ideal  shape  that  is  being  optimized,  a,  /?,  and  7  are 


angles  by  which  the  orthogonal  rotation  matrix  R  adjusts  Ipa.  Also,  A  =  RA*R  1  = 

1  T 


gives  the  final  orientation  vectors  771,  772,  and  773  that  correspond  to  the 


directions  of  a,  b,  and  c,  respectively.  Here,  ||«||F  represents  the  Frobenius  norm.  Note 
that  we  place  the  ellipsoid’s  origin  at  the  centroid  (XCl  Yc,  Zc)  of  the  voxel  object,  and  the 
density  of  all  particles  is  assumed  constant. 

The  resulting  differences  between  the  voxel  and  ideal  particles  are  considered  by  looking  at 
the  differences  in  the  volume,  surface  area,  and  inertia  tensors  after  optimizing  the  objective 
functions.  The  first  metric,  AH,  quantifies  the  difference  in  volume  (as  a  percentage)  between 
a  voxel  particle  and  its  corresponding  ellipsoid: 


0  <  AV  =  ^Vpa  Vct^  x  100  <  100  [%].  (5) 

Act  | 

The  second  metric,  AS  quantifies  the  difference  in  surface  area  (as  a  percentage)  in  the 
same  manner  as  AV  (see  Equation  (5)).  The  third  metric,  An2  =  n2  x  100  [%],  quantifies 
differences  in  orientation.  Particles  with  An2  >  5%  are  marked  in  step  vii  for  additional 
segmentation  in  step  viii  (see  orange  particle  in  Figure  4  and  description  in  subsection  II C). 

The  corresponding  idealized  pack  and  the  distribution  of  particles  for  the  voxel  pack  are 
presented  in  Figure  7,  where  e  is  a  measure  of  eccentricity  defined  as 


1 

6  “  3 


\/a*2  —  b*2  y/b*2  —  c*2  y/a2*  —  c 


*2 


+ 


+ 


(6) 


Higher  values  of  e  correspond  to  more  eccentric  particles,  and  a  perfect  sphere  corresponds  to 
e  =  0.  The  distribution  outlined  in  black  corresponds  to  the  spheres  (blue  particles),  while 
the  distribution  outlined  in  gray  corresponds  to  ellipsoids  (green  particles).  As  can  be  seen 
in  Figure  7,  the  average  eccentricity  of  the  mustard  seeds  is  small.  Thus,  they  are  referred 
to  as  spherical.  Note  that  a  substantial  reduction  of  the  dataset  size  from  3.79  x  108  voxel 
to  10,620  ellipsoids  was  achieved.  Considering  that  an  ellipsoid  is  described  by  the  centroid 
(Xc,  Yc,  Zc ),  the  semiaxes  (a,  b,  c),  and  three  angles  (a,  /3,  and  7),  we  need  9  numbers 
to  characterize  an  ellipsoidal  particle.  Thus,  95,580  numbers  are  needed  to  describe  this 
pack.  Compared  to  the  3.79  x  108  voxels,  this  dataset  is  compressed  by  approximately  four 
thousand  times. 

The  distributions  of  the  errors  in  volume,  surface  area,  and  inertia  tensor  for  all  particles 
in  the  pack  are  shown  in  Figure  8.  All  average  errors  are  below  3%.  Note  that  the  maximum 
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(a)  (b) 


FIG.  7.  Idealized  pack  and  particle  distribution  of  voxel  pack  from  Figure  6.  There  are  2276 
ellipsoidal  and  8344  spherical  particles  in  this  pack.  Compare  (a)  to  Figure  6. 

error  in  surface  area,  AS,  is  large  relative  to  AV,  because  the  weight  factor,  co,  favors 
matching  the  volume  of  the  particle  over  the  surface  area.  It  is  well  known  that  voxel 
representation  of  a  surface  is  inaccurate  and  cannot  be  improved  by  increasing  the  voxel 
resolution  [49].  Capturing  the  surface  area  of  a  voxel  object  accurately  is  nontrivial  and 
requires  fitting  a  smooth  surface  to  a  stepwise  boundary  [49,  50].  Because  a  sophisticated 
surface  representation  was  not  considered  in  this  work,  c o  was  chosen  to  be  0.8.  Nevertheless, 
only  3.99%  of  the  particles  in  this  pack  have  AS  >  5%.  In  addition,  note  that  while  n2 
is  used  to  optimize  the  orientation  of  the  idealized  shape,  this  value  is  also  a  measure  of 
how  well  the  idealized  shape  fits  the  voxel  particle  in  an  overall  sense.  For  this  pack,  only 
0.03%  of  particles  have  orientation  errors  of  AIf2  >  5%.  Note  that  using  a  traditional 
watershed-based  segmentation  algorithm  without  our  R3  mapping  (system  before  step  vii) 
would  result  in  AIf2 (max)  =  16.04%,  and  0.34%  of  particles  having  AI12  >  5%.  This  shows 
that  a  substantial  improvement  was  achieved  (AI12  (max)  was  improved  2x  and  AI12  >  5% 
was  improved  10 x). 

E.  Pack  Analysis 

Using  the  microstructure  characterization  procedures  described  above,  five  randomized 
compositions  and  one  semi-ordered  pack  of  ellipsoids  are  characterized.  These  datasets 
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FIG.  8.  (a),  (b),  and  (c)  are  the  distributions  of  particle  errors  in  volume,  surface  area,  and  inertia 
tensors  between  the  voxel  and  idealized  packs,  respectively. 


provide  a  foundational  dataset  for  analyzing  the  microstructure’s  effect  on  overall  material 
properties  that  will  be  presented  in  section  V. 


1.  Semi- ordered  Pack  of  Ellipsoids 


(a)  (b) 

FIG.  9.  Voxel  pack  of  semi-ordered  ellipsoids.  There  are  6536  ellipsoidal  particles  in  this  pack. 

The  voxel  pack  representation  after  image  processing  with  the  associated  distribution  of 
particles  is  presented  in  Figure  9,  where  the  size  of  the  cylindrical  domain  analyzed  has  a 
diameter  of  D  —  59.44  mm  and  a  height  of  H  =  47.16  mm.  This  pack  had  a  particle  volume 
fraction  of  cp’B  =  0.6841  before  particle  segmentation  and  was  reduced  to  cp,A  =  0.6655  after 
image  segmentation  steps,  corresponding  to  a  volume  loss  of  £jp  =  2.71%.  The  VOI  for  this 
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sample  is  -Dyoi  =  45.24  mm  and  //yoi  =  33.03  mm.  When  modeling  the  voxel  particles 
as  idealized  ellipsoids,  the  mean  errors  of  AV,  AS  and  An2  are  0.70%,  4.48%,  and  2.20%, 
respectively.  The  maximum  errors  in  AV,  AS  and  AII2  are  2.69%,  20.98%,  and  11.65%, 
and  24.91%  and  0.69%  of  all  particles  have  AS  >  5%  and  An2  >  5%. 

In  Figure  9(a),  it  can  be  seen  that  the  particles  are  aligned  in  the  y-direction.  The 
orientation  of  an  ellipsoidal  particle  is  described  by  two  angles  with  respect  to  the  basis 
vector  of  the  longest  semi-axis  f]i :  the  azimuthal  angle  (projected  orientation  onto  the  xy 
plane)  and  the  altitude  angle  (with  respect  to  the  z  axis).  These  angles  are  referred  to  as 
9a  and  (pa  (see  Figure  10)  and  are  found  from  the  orientation  of  an  idealized  shape.  The 
distributions  of  the  azimuthal  and  altitude  angles  for  the  semi-ordered  pack  is  presented  in 
Figure  11.  The  distribution  of  cpa  shows  that  the  ellipsoids  are  more  likely  to  lie  horizontally 
(<pa  =  0°)  than  in  the  unstable  vertical  position  (<pa  =  90°),  while  the  distribution  of  9a 
shows  that  the  ellipsoids  are  preferably  aligned  in  the  y-direction  (9a  =  —90°  =  90°),  which 
can  also  be  observed  in  Figure  9(a).  This  anisotropic  morphology  will  have  a  critical  impact 
on  the  overall  thermal  conductivity  estimates  presented  in  section  V. 


Z 


FIG.  10.  Definition  of  ellipsoid’s  orientation  with  respect  to  the  largest  semiaxis. 


2.  Randomized  Mixtures 

Five  compositions  at  0%,  25%,  50%,  75%  and  100%  weight  of  spheres  to  ellipsoids  are 
considered.  These  five  compositions  are  created  by  weighing  the  appropriate  amount  of  each 
constituent  and  are  then  packed  and  tomographically  characterized  two  times  to  confirm  the 
repeatability  of  the  steps  described  in  this  section  (We  refer  to  these  two  scans  as  higher 
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FIG.  11.  Distribution  of  particle  orientations  for  semi-ordered  pack  shown  in  Figure  9(a). 


and  lower  volume  fraction  packs  depending  on  their  packing  density).  Table  I  shows  infor¬ 
mation  about  each  voxel  pack  considered  including  the  weight  of  each  constituent,  sizes  of 
the  datasets,  and  the  error  introduced  due  to  the  image  segmentation  steps.  The  largest 
volume  losses,  e  j y> ,  occur  for  the  spherical  particle  packs.  However,  all  volume  losses  are 
less  than  3.66%.  Note  that  the  height  of  the  volume  analyzed,  H ,  increases  with  increasing 
weight  %  of  spherical  particles,  because  the  spherical  particles  are  less  dense  than  the  el¬ 
lipsoidal  ones.  This  can  be  also  observed  in  Figure  2,  since  darker  gray  corresponds  to  less 
dense  material.  The  resulting  volume  fractions  of  all  randomized  compositions  considered 
in  this  work  before  and  after  image  segmentation,  cvp’B  and  cp’A,  are  presented  in  Figure 
12(a).  Figure  12(b)  shows  the  volume  fraction  of  each  constituent  for  the  corresponding 
composition.  In  these  figures,  the  two  samples  for  each  composition  are  represented  by  the 
two  data  points,  and  the  lines  are  the  means.  All  differences  in  the  total  particle  volume 
fractions  after  segmentation,  between  the  two  packs  of  each  composition  analyzed  (higher 
and  lower  volume  fraction  packs),  are  below  0.008  (maximum  difference  is  for  25%  mixture, 
see  Figure  12(a)).  The  maximum  differences  between  the  phase  volume  fractions,  cve  and 
Cg,  are  0.0191  and  0.0156,  respectively,  for  the  50%  mixture  (see  Figure  12(b)).  Recall  that 
cp,A  =  Cg  +  cvs.  As  the  particle  volume  fractions,  cp,B ,  can  theoretically  be  in  the  range  [0,1], 
these  differences  between  the  two  data  points  at  each  composition  are  small.  The  number  of 
particles  per  volume  is  plotted  for  each  of  the  five  packs  in  Figure  13,  where  Np  is  the  number 
of  particles  in  the  VOI.  The  number  of  total  particles  increases  with  increasing  weight  %  of 
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Pack  [%  wt. 

I  0 

25 

50 

75 

100 

we 

[g] 

100.001  74.998 

50.003 

25.007 

0.000 

ws 

[g] 

0.000 

24.997 

50.003 

75.002 

100.000 

D  [mm]  - 57.67 


Dvoi  [mm] -  46.57  - 

Higher  Volume  Fraction  Packs 


H  [mm] 

41.43 

42.74 

48.51 

50.04 

51.42 

Hvoi  [mm 

]  30.33 

31.64 

37.41 

38.93 

40.32 

SIP  [%] 

2.23 

2.85 

2.60 

2.72 

3.66 

Lower  Volume  Fraction  Packs 

H  [mm] 

41.22 

43.10 

46.22 

49.62 

49.76 

Hyoi  [mm 

]  30.12 

31.99 

35.12 

38.52 

38.66 

SIP  [%] 

2.22 

2.72 

2.46 

2.66 

3.52 

TABLE  I.  Description  of  randomized  packs.  Ws  and  We  are  the  weights  of  each  particulate  phase 
(s  subscript  for  sphere  and  e  subscript  for  ellipsoid). 


FIG.  12.  Volume  fraction  plots  for  various  compositions  of  the  surrogate  system.  Note  that  two 
samples  for  each  composition  were  studied  and  are  referred  to  as  higher  volume  fraction  packs  (top 
points)  and  lower  volume  fraction  packs  (bottom  points). 

spherical  particles,  because  the  volume  of  a  spherical  particle  is  smaller  than  the  volume  of 
an  ellipsoidal  particle  as  evident  in  the  distribution  of  particles  for  the  50%  pack  in  Figure 
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FIG.  13.  Number  of  particles  per  cm3. 

7(b).  The  mean  and  max  errors  between  the  voxel  particles  and  their  corresponding  ideal 
ellipsoidal  representation  are  presented  in  Table  II.  All  mean  errors  for  each  pack  are  below 
5%.  As  discussed  earlier,  the  noticeable  maximum  errors  in  the  surface  area  are  related  to 
errors  in  representing  the  surface  of  voxel  particles.  Moreover,  the  noticeable  maximum  er¬ 
rors  for  II2  exist  for  a  small  number  of  particles  that  do  not  match  an  ellipsoid  well.  Overall, 
the  low  errors  achieved  throughout  the  microstructure  characterization  process  show  that 
the  sample  preparation,  data  acquisition,  and  image  processing  steps  produce  scientifically 
sound  and  repeatable  results  with  well  controlled  errors. 

The  mixture  of  randomized  ellipsoids  (0%  wt.  spherical  particles)  has  a  volume  fraction 
consistent  with  packing  simulations  of  monodisperse  ellipsoids  reaching  a  maximally  random 
jammed  state  [51].  For  the  average  size  ellipsoidal  particle  in  this  work  (a  =  3.04  mm, 
b  =  1.11  mm,  and  c  =  0.83  mm),  a  maximum  packing  volume  fraction  can  be  extrapolated 
to  be  0.68  from  the  work  of  Chaikin  et  al.  [51]  ( cvp'B  =  0.66  in  this  work).  It  is  expected 
that  this  volume  fraction  should  be  similar  as  the  ellipsoidal  particles  have  a  small  range  in 
size  as  seen  in  Figure  7.  The  volume  fraction  of  maximally  jammed  monodisperse  spheres 
in  the  work  of  Chaikin  et  al.  [51]  is  0.64,  and  the  volume  fractions  for  the  compositions 
with  100%  spheres  considered  in  this  work  is  higher  (cp,B  =  0.70  from  Figure  12)  due  to  the 
polydispersity  as  evidenced  by  the  range  in  sizes  in  Figure  7.  Note  that  this  volume  fraction 
falls  below  0.74,  the  assumed  maximum  achievable  volume  fraction  of  monodisperse  spheres 
first  proposed  by  Keplar  and  supported  by  many  thereafter  [8] .  When  considering  all  of  the 
randomized  compositions,  there  is  a  monotonically  increasing  trend  in  volume  fraction,  cvp. 
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Pack  [%  wt.] 
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25  50 

75 

100 

Higher  Volume  Fraction  Packs 

AV  (mean) 

0.65 

0.50  0.41 

0.33 

0.41 

AV  (max) 

3.63 

3.29  3.29 

3.27 

3.69 

AS  (mean) 

4.11 

3.16  2.56 

2.03 

2.62 

AS  (max) 

34.22  28.55  28.62  28.15  35.28 

An2  (mean)  1.09  0.56  0.38  0.22  0.21 
An2  (max)  14.17  8.41  8.37  5.68  18.69 


Lower  Volume  Fraction  Packs 


AV  (mean) 

0.63 

0.52 

0.35 

0.32 

0.40 

AV  (max) 

3.69 

3.52 

3.50 

3.22 

3.54 

AS  (mean) 

3.99 

3.31 

2.17 

2.04 

2.50 

AS  (max) 

35.81  32.26  26.05  27.56  32.71 

An2  (mean) 

1.05 

0.57 

0.36 

0.22 

0.21 

An2  (max) 

6.47 

12.79 

11.43  4.97 

10.09 

TABLE  II.  Summary  of  errors  between  voxel  particles  and  idealized  ellipsoids.  The  unit  of  all 
values  is  %.  Less  than  0.1%  of  all  particles  considered  in  these  mixtures  have  AII2  >  5%,  and  less 
than  14.3%  of  all  particles  have  AS  >  5% 


It  is  known  that  when  the  scales  of  two  inclusions  are  significantly  different,  the  mixtures 
will  have  higher  volume  fractions  than  the  maximum  volume  fractions  of  a  single  constituent 
for  monomodal  spheres  [52],  However,  as  the  smallest  semiaxes  of  the  ellipsoids  are  similar 
to  the  average  radius  of  the  spherical  particles,  no  such  effect  is  observed. 

The  distributions  of  particle  orientations  for  the  randomized  pack  with  100%  weight  of 
ellipsoids  is  shown  in  Figure  14.  Notice  that  there  is  no  preference  of  the  randomized  pack’s 
ellipsoids  azimuthal  angle,  6a  (compare  with  Figure  11  for  semi-ordered  pack).  Moreover, 
there  is  a  low  probability  that  particles  are  suspended  in  the  unstable  vertical  position, 
( f>a  =  0°.  This  randomness  in  the  azimuthal  angle,  9a,  confirms  the  randomizing  effects  in 
the  sample  preparation  steps.  The  ellipsoidal  particles  in  the  other  randomized  packs  exhibit 
similar  trends  in  their  orientations  and  are  not  reported  here. 
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0.025 


FIG.  14.  Distribution  of  orientations  for  ellipsoidal  particles  in  the  pack  with  0%  weight  of  spherical 
particles. 


When  comparing  the  semi-ordered  pack  of  ellipsoids  to  the  randomized  one,  the  volume 
fraction  of  the  semi-ordered  ellipsoid  pack  is  higher.  This  is  in  agreement  with  studies  of 
ordered  packs  of  monodisperse  ellipsoids,  as  they  are  known  to  have  a  maximum  particle 
volume  fraction  of  0.74  —  0.77  [53].  When  comparing  the  altitude  angle,  (f>a,  distributions 
between  the  semi-ordered  and  randomized  pack,  there  is  a  higher  probability  of  the  semi- 
ordered  pack  particles  existing  in  the  xy  plane  (</>a  =  0°)  and  a  preferred  direction  for  6a 
(see  Figures  11  and  14). 

Let  us  note  that  the  microstructure  characterization  method  presented  in  this  work  can 
be  used  in  a  variety  of  applications  as  well  as  in  validation  of  packing  simulations  [10,  54].  In 
this  work,  the  datasets  provide  a  foundational  set  for  considering  the  microstructure’s  effect 
on  overall  properties,  and  the  effect  of  these  packing  configurations  on  the  overall  anisotropic 
thermal  conductivity  tensor  will  be  explored  in  the  remainder  of  this  paper. 


III.  THEORETICAL  OVERVIEW  OF  OVERALL  ANISOTROPIC  PROPERTIES 

Using  the  Hashin-Shtrikman  variational  principle  and  assuming  that  particulate  com¬ 
posites  are  homogeneous  and  ergodic,  Willis  formulated  a  second-order  approximation  of 
bounds  on  elastic  and  conductivity  constants  in  which  the  two-point  probability  functions 
are  embedded  within  the  mathematical  formulation  directly  [28].  In  Willis’  work,  assump¬ 
tions  about  the  configuration  of  the  particles  were  made  when  computing  overall  properties. 
In  this  work,  upper  and  lower  anisotropic  bounds  as  well  as  anisotropic  self-consistent  esti- 
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mates  are  computed  for  the  first  time  to  the  best  of  our  knowledge  without  assumptions  of 
statistical  isotropy  nor  ellipsoidal  or  any  other  material  symmetry  for  complex  tomograph- 
ically  characterized  polydisperse  microstructures.  A  simple  overview  of  theory  is  discussed 
next. 


A.  n-point  Probability  Functions 


As  described  in  [16,  42],  the  n-point  probability  functions  are  derived  using  a  phase 
indicator  function  at  a  position  x  in  a  sample  a  of  an  ensemble  space  £ : 

{1  if  x  in  phase  r 

(7) 

0  otherwise. 

An  ensemble  is  a  collection  of  material  samples  being  considered.  The  ensemble  average  is 
given  by 

Xr(x)  =  J  Xr(x',  ot)p(a)da,  (8) 

where  p(a)  is  the  probability  density  function  of  a  in  £.  The  n-point  probability  function, 
Srir2... rn{xi,  x2,  ■  ■  ■  ,  xn) ;  is  defined  as 

Srs---q(xli  x2i  j  xn)  Xr(x  1  )Xs(x 2)  '  '  '  Xg(Xn)-  (9) 


It  represents  the  probability  of  phases  r,  s,  ■  ■  ■  ,  q  existing  at  points  xi,  x2,  •  •  •  ,  xn,  simul¬ 
taneously.  In  general,  the  probability  functions  for  a  heterogeneous  material  are  spatially 
complex.  I11  this  work,  a  second-order  model  for  computing  bounds  on  overall  anisotropic 
thermal  conductivity  is  considered. 

For  statistical  homogeneous  materials,  where  the  probability  functions  are  translationally 
invariant,  the  one-point  probability  function  attains  a  constant  value,  Sr(x)  =  c,  and  the 
two-point  function  simplifies  to  Srs(x,  x')  =  Srs(x  —  x').  For  statistically  homogeneous  sys¬ 
tems,  it  is  meaningful  to  define  volume  averages.  When  assuming  ergodicity  of  homogeneous 
systems,  ensemble  averaging  is  equivalent  to  volume  averaging  in  the  infinite  volume  limit 


Srs---q(x  1;  x2i  '  '  '  1  xn) 


lim  7T  /  XrOl  -  l)Xs{x 2  -  0  ■  ■  ■  Xq(xn  ~  l )  dfl, 

il— >00  i  l  / o 


(10) 
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where  l  is  a  translation  vector,  and  is  the  volume  of  the  domain.  For  homogeneous  and 
ergodic  systems,  the  one-point  probability  function  is  equivalent  to  the  volume  fraction  of 
phase  r  (Sr(x)  =  cr).  The  two-point  probability  functions  of  a  statistically  homogeneous 
and  isotropic  system  (rotationally  and  translationally  invariant)  are  defined  as  Srs(x,x')  = 
SVs(|zi  —  £c'|),  where  the  function  only  depends  on  the  distance  between  two  points.  For 
homogeneous  systems  with  no  long-range  order,  we  observe  two  limit  cases  in  the  pointwise 
sense  of  the  two-point  probability  functions,  which  can  be  expressed  as 

{cr8rs  as  x  —  x’  — y  0 

(11) 

crcs  as  x  —  x'  — >•  oo, 

where  8rs  is  the  Kronecker  delta.  More  information  on  the  n-point  probability  functions  can 
be  found  in  work  by  Torquato  [55]  and  Beran  [56].  In  this  work,  statistically  anisotropic 
systems  are  considered,  and  n-point  probability  functions  are  used  when  addressing  the 
issues  of  material  order. 

B.  Hashin-Shtrikman  Variational  Principle 

The  Hashin-Shtrikman  variational  principle  is  employed  to  obtain  bounds  on  thermal- 
mechanical  properties  of  an  anisotropic  material  as  described  in  [27,  28].  This  paper  focuses 
on  computing  the  overall  bounds  and  self-consistent  estimates  for  the  anisotropic  thermal 
conductivity  tensor.  Note  that  this  formulation  is  similar  to  theory  for  overall  linear  elastic 
constants  [42], 

It  is  assumed  that  the  material  follows  a  constitutive  relation  of  q  =  d w(Q)/dQ,  where 
w(Q)  is  an  energy  density  function,  q  is  the  heat  flux  vector,  and  Q  =  —  VT  is  the  negative 
temperature  gradient  for  a  temperature  T  at  a  point  x.  When  the  material  obeys  Fourier’s 
law,  q  =  k-Q  (linear  constitutive  relation),  the  energy  density  function  is  defined  as  w(Q)  = 
(Q  •  k  ■  Q)/2  >  0,  where  n  is  the  symmetric  positive  semi-definite  second  order  thermal 
conductivity  tensor.  This  symmetric  positive  semi-definiteness  requirement  is  in  agreement 
with  Onsager’s  reciprocity  theorem  [57] .  Thus,  the  generally  anisotropic  conductivity  tensor 
computed  in  this  work  has  six  independent  components.  A  complementary  energy  density, 
w*(q),  can  also  be  defined,  where  w*(q)+w(Q)  =  q  -Q  and  Q  =  dw*(q)/dq.  For  a  material 
following  Fourier’s  law,  w*  —  (q  ■  p  ■  q)/ 2  >  0,  where  p  represents  the  symmetric  positive 
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semi-definite  second  order  thermal  resistivity  tensor.  Note  that  p  =  n~l.  The  internal 
energy  and  complementary  internal  energy  of  the  system  are  defined  as  E  =  fQw  df2  and 
E*  =  fQ  w*  df2,  respectively. 


T  =  T  on  an 


To  =  T  on  an 


t'  =  o  on  an 


qn  —  qn  on  an 

General  Body 


q°  =  qn  on  <9f 1  q'n  =  0  on  dQ 

Comparison  Medium  Polarization  Field 


FIG.  15.  Body  decomposition  with  prescribed  Dirischlet  and  Neumann  boundary  conditions,  re¬ 
spectively,  corresponding  to  the  lower  bound  and  upper  bound  formulations. 


Hashin  and  Shtrikman  [27]  proposed  the  body  decomposition,  as  shown  in  Figure  15. 
The  top  half  of  this  figure  shows  the  body  decomposition  with  a  prescribed  temperature 
held  on  the  boundary,  T,  which  is  used  to  compute  the  lower  bound  (LB)  of  the  overall 
77.  The  bottom  half  of  this  figure  shows  the  decomposition  for  a  body  with  a  prescribed 
normal  heat  flux,  qn,  and  this  formulation  is  used  for  computing  the  lower  bound  of  the 
overall  p,  corresponding  to  the  upper  bound  (UB)  of  77.  Note  that  quantities  with  •  denote 
an  overall  or  volume  averaged  quantity,  i  —  1/Q  f  •  dfh  The  governing  equations  of  the 
corresponding  elliptic  boundary  value  problems  are  given  by 


LB:  UB: 


V  ■  [k0  ■  Q  +  p\  =0 

in  12, 

<1 

o 

in 

12, 

(13) 

(12) 

p  —  [k  —  k0]  •  Q  =  0 

in  fl, 

t  -  [p  *  po]  ■  q  =o 

in 

12, 

T'  =  0 

cf 

Pl 

o 

q  ■  n  =0 

on 

<912, 

where  k,q  is  the  thermal  conductivity  tensor  of  a 

comparison  medium, 

and 

p  is 

the  heat 

flux  polarization  tensor.  To  and  T'  represent  the  homogeneous  and  fluctuation  temperature 
fields,  respectively,  q'  is  the  heat  flux  fluctuation  held,  po  is  the  comparison  medium  resis¬ 
tivity  tensor,  t  is  the  temperature  gradient  polarization  held,  qn  =  q  ■  n  is  the  normal  heat 
hux,  and  n  is  the  unit  normal  vector.  Formulations  equivalent  to  the  ones  described  by  the 


23 


strong  forms  (12)  and  (13)  can  be  obtained  by  minimizing  the  functionals 


LB: 


2 Jr{p)=  I  p-[k-k0\  1  •  p  +  p  ■  I  T(x,  x ')  •  [p(aj')  -  p]dfix/  -2 p  Qc 
J  Q  _  J  Q  / 


dfi, 


(14) 


and 


UB: 


2E*(t)  =  /  t  •  [p  —  p0]  1  •  t  +  t  •  /  [/(#,  a;')  •  [£(&')  —  £]d£)x/  —  2t  ■  q0 
A)  _  «/ r2„/ 


dfi, 


(15) 


for  the  lower  and  upper  bound  formulations,  respectively.  r(ic,  x')  is  the  second  order  linear 
operator  related  to  the  Green’s  function  solution.  Q0  denotes  the  negative  temperature 
gradient  in  a  homogeneous  comparison  medium,  and  qo  is  the  homogeneous  heat  flux.  Note 
that  the  energy  functional  is  defined  as  E  —  Eq  —  E,  where  E0  is  the  internal  energy  of 
the  comparison  medium.  T*  =  E*  —  E*,  where  Eq  is  the  complementary  internal  energy 
of  the  comparison  medium.  For  the  upper  bound  formulation,  U  is  the  second  order  linear 
operator  related  to  the  Green’s  function  solution  (U  =  k,  —  k  ■  T  •  k). 

To  find  second-order  estimates  to  the  stationary  point  of  Equations  (14)  and  (15),  trial 
fields  of  p  and  t  are  assumed  piecewise-constant  in  each  phase  r  ( p *  =  J2r= o  XrPr  and 
t*  =  Ylr=oXrtr  f°r  a  system  of  n  particulate  phases  and  a  matrix  corresponding  to  r  =  0). 
The  following  systems  of  equations  are  obtained  after  discretization  of  the  polarizations 
fields,  p  and  t,  averaging  over  the  ensemble  space  £  and  using  the  calculus  of  variations 


LB: 


n  „ 

cr(nr  —  Ko)^1  ■  Qr  +  /  r°°(a3  —  x')(Srs(x  —  x’)  —  crcs)  ■  qsd£lx' 

s=0 


CrQi 


(16) 


and 


UB: 


n  « 

Cr(pr  ~  Pc)'1  •  tr  +  /  U°°(x  —  x')(Srs(x  —  x')  —  CrCs)  ■  ts  dnx, 

s=0  -’nx' 


—crq. 


(17) 


Here  cr  and  cs  are  the  volume  fractions  of  phase  r  and  s,  respectively,  Kr  is  the  thermal 
conductivity  tensor  of  phase  r,  while  Q  is  the  overall  negative  temperature  gradient.  r°° 
and  U°°  are  related  to  the  Green’s  function  solution  for  an  infinite  body.  For  the  upper 
bound  formulation,  p.r  is  the  resistivity  tensor  of  phase  r,  while  q  is  the  overall  heat  flux. 
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These  systems  of  equations  are  solved  for  the  piecewise-constant  fields,  pr  and  tr,  to 
find  the  mean  polarization  fields,  p  =  Y2r=ocrPr  and  ^  =  J2r=oc^r-  The  volume  averaged 
constitutive  relations, 


LB:  UB: 

(18) 

q  =  KQ  =  K0Q  +  p,  Q  =  pq  =  p0q  +  t, 

are  used  to  calculate  the  overall  77  and  p,  respectively.  When  Kr  —  re0  is  the  smallest  positive 
semi-definite  matrix  for  all  phases  r  in  the  lower  bound  formulation,  the  resulting  77  is  the 
second-order  lower  bound  of  the  conductivity  tensor.  Meanwhile,  when  pr  —  po  is  the  smallest 
positive  semi-definite  matrix  for  all  phases  r  in  the  upper  bound  formulation,  the  second- 
order  lower  bound  of  the  resistivity  tensor  is  obtained  (upper  bound  of  the  conductivity 
tensor,  77  =  p-1). 

Another  widely  used  model  for  computing  overall  material  properties  is  the  self-consistent 
estimate  (SC).  As  described  in  Willis’  work  [28],  self-consistent  estimates  are  calculated 
by  minimizing  the  volume  averaged  energy  and  complementary  energy  functionals  \J-\  = 
\E0  —  E\  and  \E  |  =  \E^  —  E*\.  This  leads  to  E  =  E0  (k  =  n0)  and  E*  =  E$  (p  —  Po)- 
With  this  in  mind,  the  following  objective  functions  are  minimized  in  order  to  compute 
self-consistent  estimates 


SC-L: 


WL  — 


SC-U: 


r\u  — 

li5C  ~ 


[Po 


(19) 


The  integrals  presented  in  Eq.  (16)  and  (17)  and  consequently  in  (19),  which  are  products 
of  the  tensors  related  to  the  Green’s  function  solution,  T00  and  U°°,  and  the  second-order 
probability  function,  are  strenuous  especially  near  the  origin  due  to  the  singularity  of  T00 
and  U°°.  Note  that  for  statistically  isotropic  systems  and  systems  with  aligned  ellipsoids, 
these  integrals  can  be  calculated  analytically  leading  to  a  closed  form  solution  of  the  overall 
conductivity  tensor.  Closed  form  solutions  using  these  assumptions  have  been  presented  in 
work  by  Willis  [28]  and  Weng  [58],  just  to  name  a  few.  This  work  is  not  limited  by  these 
assumptions,  and  the  integral  kernels  are  calculated  numerically  using  the  adaptive  sparse- 
grid  Smolyak  integration  method  with  hierarchical  basis  that  was  developed  in  our  previous 
work  [42], 
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IV.  STATISTICAL  CHARACTERIZATION 


After  the  voxel  data  is  processed,  the  packs  are  statistically  characterized  by  computing 
one-  and  two-point  probability  functions  directly  from  the  voxel  data  using  our  parallel 
statistical  sampling  code,  STAT3D  [16].  Note  that  sampling  in  M3  on  idealized  packs  as 
done  in  [16]  is  also  possible  and  would  result  in  nearly  identical  results  as  shown  in  [22], 
We  do  not  present  this  verification  here.  This  section  presents  the  statistical  functions 
for  compositions  described  in  section  HE.  The  semi-ordered  pack  and  one  realization  of 
each  randomized  composition  (higher  volume  fraction  packs  as  shown  in  Figure  12(a))  are 
statistically  characterized. 

As  presented  in  [16],  a  spherical  sampling  template  is  used  to  compute  n-point  probability 
functions.  After  a  convergence  study,  the  sampling  template  was  discretized  with  2000  radial 
points  and  32  circumferential  points  and  was  randomly  thrown  into  the  sample  10  million 
times.  Note  that  this  very  high  sampling  frequency  is  important  for  resolving  the  probability 
functions.  In  order  to  quantify  the  errors  in  the  statistical  functions,  we  consider  the  limit 
cases  of  the  two-point  probability  functions  as  given  by  Eq.  (11).  The  following  error 
measures  at  the  origin  and  at  a  far  distance  are  used 


(20) 


(21) 


where  std(«)  is  the  standard  deviation. 

If  a  pack  is  perfectly  isotropic,  the  standard  deviations  in  the  two-point  probability  func¬ 
tions  at  \x  —  x'\  are  zero  (std(SVs(cc  —  a/)l|!E-a:,|)  =  0).  For  random  statistically  homogeneous 
mixtures,  short  range  statistical  anisotropy  can  exist  due  to  both  the  shape  of  the  particles 
and/or  the  configuration  of  the  particles,  and  no  long  range  order  is  expected  for  random¬ 
ized  systems.  A  specific  case  of  statistically  anisotropic  systems,  considered  in  previous 
works  [27,  28],  is  when  the  two-point  probability  functions  exhibit  ellipsoidal  symmetry, 
which  is  known  to  occur  for  ordered  packs  of  randomly  positioned  and  uniformly  aligned 
monodisperse  ellipsoids.  Ellipsoidal  symmetry  means  that  the  function  values  of  Srs(x  —  x') 
are  constant  on  contours  of  an  ellipsoid, 


(aj  —  x’Y  B  1(x  —  x1), 


(22) 
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where  B  is  a  positive  definite  matrix  whose  eigenvalues  are  the  squares  of  an  ellipsoid’s 
semiaxes.  It  will  be  shown  that  the  statistical  descriptors  of  the  polydisperse  packs  consid¬ 
ered  in  this  work  do  not  exhibit  ellipsoidal  symmetry.  When  ellipsoidal  symmetry  does  not 
exist,  closed-form  solutions  of  the  overall  properties  as  presented  in  [28]  and  many  papers 
thereafter  can  not  be  obtained. 

A.  Semi-ordered  Pack 


(a) 


FIG.  16.  Mean  and  standard  deviation  of  two-point  probability  functions  for  semi-ordered  pack  of 
ellipsoids  from  Figure  9. 


In  what  follows,  the  subscripts  m,  s  and  e  are  used  to  refer  to  the  matrix,  spherical 
particles  (mustard  seeds),  and  ellipsoidal  particles  (rice  grains),  respectively.  The  one-point 
probability  functions  for  the  semi-ordered  mixture  of  ellipsoids  (see  Figure  9)  are  cm  =  0.3345 
and  ce  =  0.6655.  The  mean  and  standard  deviations  of  the  two-point  probability  functions 
are  presented  in  Figure  16.  The  maximum  errors,  as  defined  by  Equations  (20)  and  (21) 
for  all  two-point  probability  functions  of  this  mixture  are  £srs(o)  —  0.0059%  and  £sra( oo)  = 
0.304%,  respectively.  As  can  be  seen  in  Figure  16(a),  the  mean  two-point  probability  values 
reach  local  extrema  near  the  average  semiaxis,  (a  +  b  +  c)/ 3  =  1.66  mm  (vertical  gray  dotted 
line)  and  the  largest  semiaxis  (vertical  black  dashed  line).  Note  that  the  average  semiaxes  of 
ellipsoids  in  this  work  are  a  =  3.04  mm,  b  =  1.111  mm,  and  c  =  0.83  mm.  Considering  the 
mean  function  values  and  the  standard  deviations  as  ( x  —  x ')  oo  (Figure  16),  the  functions 
are  considered  converged  when  \x  —  x'\  =  8  mm,  resulting  in  a  characteristic  material  length 
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scale  of  16  mm  (twice  the  radius) .  The  characteristic  material  length  scale  can  be  considered 
as  the  optimal  size  of  a  Representative  Unit  Cell  (RUC)  [16].  This  length  scale  is  ~  20  x 
the  smallest  semi-axis  and  ~  5x  the  largest  semi-axis.  The  characteristic  material  length 
scale  is  important  in  this  work,  as  the  saturation  of  the  probability  functions  guarantee  the 
convergence  of  the  integrals  in  Equations  (16),  (17)  and  (19)  as  ( x  —  x')  — >■  oo  [28].  Also  of 
interest  in  Figure  16(b)  is  where  the  maximum  standard  deviations  occur,  since  the  std(S'rs) 
represents  how  anisotropic  the  system  is.  The  largest  peak  occurs  near  a  distance  from  the 
origin  of  the  smallest  semiaxis  (c  =  0.83  mm),  the  second  peak  near  the  average  semiaxis 
(dotted  gray  line),  and  the  deviations  diminish  after  the  distance  near  the  largest  semiaxis 
(dashed  black  line).  Note  that  these  deviations  are  related  to  a  combination  of  the  particle 
shape,  the  polydispersity,  and  the  configuration.  While  these  deviations  are  less  than  6%  of 
the  function  values,  this  statistical  anisotropy  will  lead  to  larger  macroscopic  anisotropy  of 
the  thermal  conductivity  tensor,  because  the  anisotropy  in  morphology  is  amplified  by  the 
local  anisotropic  material  constants. 

One  of  the  two-point  probability  functions,  See,  is  visualized  in  Figure  17.  Notice  again 
that  local  minima  exist  near  distances  from  the  origin  to  the  average  semiaxis  and  the 
largest  semiaxis  lengths.  The  anisotropic  nature  of  the  two-point  probability  function  in  the 
xy  and  yz  plane  (Figure  17(c)  and  17(d),  respectively)  is  due  to  the  preferred  alignment  of 
the  particles  in  the  y  direction  (see  Figure  11).  While  the  general  shape  of  this  statistical 
function  is  ellipsoidal,  notice  that  S 22  does  not  exhibit  complete  ellipsoidal  symmetry  (given 
by  Equation  (22)),  thus  prohibiting  the  closed  form  integration  given  in  [28].  Near  spherical 
symmetry  occurs  only  in  the  xz  plane  (Figure  17(b))  since  the  smallest  semiaxes  of  the 
averaged  size  ellipsoid  are  similar  ( b  =  1.11  mm,  and  c  =  0.83  mm).  Note  that  for  a  semi- 
ordered  pack  (based  on  our  packing  method  meaning  not  randomly  positioned  and  uniformly 
aligned  inclusions)  of  ellipsoids  in  which  all  semi-axes  are  significantly  different,  ellipsoidal 
symmetry  would  not  be  present  in  any  plane. 


B .  Randomized  Packs 

The  mean  and  standard  deviations  of  the  two-point  probability  functions  for  selected  ran¬ 
domized  compositions  are  considered  in  Figures  18(a)-(f).  The  maximum  errors,  as  defined 
by  Equations  (20)  and  (21),  for  all  two-point  probability  functions  of  the  five  randomized 


x  [mm]  y  [mm] 

(c)  (d) 


FIG.  17.  (a)  3D  representation  of  See(x  —  x').  (b)-(d)  show  that  See(x  —  x')  in  the  xz,  xy  and  yz 
planes,  respectively.  Note  that  the  scale  of  the  colorbar  is  adjusted  to  highlight  where  local  minima 
of  this  function  exist. 

compositions  are  £srs{o)  =  0.104%  and  £srs( oo)  =  4.146%,  respectively.  Note  that  the  spike 
in  the  standard  deviations  at  the  origin  in  Figure  18(f)  (less  than  1%  of  the  function  value) 
occurs  at  the  voxel  resolution  of  the  dataset,  \x  —  x'\  k,  0.07  mm.  This  is  a  numerical 
artifact  due  to  the  discontinuities  associated  with  the  voxelization  of  data.  In  Figure  18(d), 
the  standard  deviations  for  the  50%  pack  as  ( x  —  x')  — >■  oo  are  not  fully  saturating  for 
some  of  functions.  This  indicates  that  larger  samples  might  be  needed  for  better  resolution. 
Nevertheless,  these  functions  are  not  significantly  misrepresenting  the  statistical  anisotropy 
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in  the  system,  since  the  standard  deviations  as  (x  —  x')  — >■  oo  are  on  the  same  order  as 
the  standard  deviation  peaks  near  the  origin.  For  the  25%  and  75%  compositions,  the  stan¬ 
dard  deviations  as  (a;  —  x')  — >  oo  are  similar  to  those  for  the  50%  mixture.  The  general 
trend  shows  that  the  maximum  standard  deviations  for  the  two-point  probability  functions 
of  the  mixtures  decrease  as  the  amount  of  spheres  in  the  system  increases  (Figures  18(b), 
(e),  and  (f)),  indicating  that  the  spherical  particles  tend  to  organize  randomly  (isotropi¬ 
cally).  For  the  100%  mixture  of  spheres,  the  deviations  are  essentially  zero  (Figure  18(f)), 
indicating  the  pack  is  nearly  isotropic.  This  trend  in  statistical  anisotropy  will  be  reflected 
in  the  macroscopic  anisotropy  of  the  thermal  conductivity  tensor.  Once  more,  individual 
local  minima  and  corresponding  high  standard  deviations  are  due  to  particle  shape  and  the 
spatial  organization. 


Considering  both  the  deviations  as  (a;  —  x')  — y  oo  and  where  the  two-point  probability 
functions  saturate,  characteristic  material  length  scales  are  determined  for  the  mixtures.  The 
mixtures  with  only  one  particulate  phase  (0%  and  100%  weight  of  spherical  particles)  have  a 
characteristic  material  length  scale  around  16  mm,  as  the  functions  converge  at  \x  —  x'\  =  8 
mm.  This  same  length  scale  is  attributed  to  the  similarity  in  the  average  semiaxis  length 
of  the  ellipsoids,  (a  +  b  +  c)/3  =  1.66  mm,  and  the  radius  of  the  average  sphere,  0.94  mm. 
The  three  phase  mixtures  (25%,  50%  and  75%  weight  of  spherical  particle  packs)  have  larger 
characteristic  material  length  scales,  20  mm.  This  larger  length  scale  is  expected  because 
more  interactions  must  be  captured  as  polydispersity  grows. 


When  comparing  the  probability  functions  of  the  randomized  pack  of  100%  ellipsoids 
to  the  semi-ordered  pack  of  ellipsoids,  the  amount  of  statistical  anisotropy  reflects  the  dif¬ 
ferences  in  configuration.  The  maximum  standard  deviation  for  the  randomized  pack  of 
ellipsoids  is  std(S'rs)  =  0.0049  (see  Figure  18(b)),  which  is  smaller  than  the  maximum  stan¬ 
dard  deviation  of  the  semi-ordered  pack  std(S'rs)  =  0.0072  (Figure  16(b)),  and  indicates  less 
statistical  anisotropy.  Also  note  that  there  are  only  two  maxima  in  the  standard  deviations 
as  compared  to  three  for  the  semi-ordered  pack  (compare  Figures  16(b)  and  18(b)).  As  there 
is  no  apparent  alignment  in  the  xy  plane  for  the  randomized  pack  (Figure  14),  the  two-point 
probability  functions  reflect  this  behavior  and  smooth  out  these  standard  deviation  peaks. 
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(e) 


(f) 


FIG.  18.  Mean  and  standard  deviations  of  two-point  probability  functions  for  randomized  packs 
of  0%  ((a)-(b)),  50%  ((c)-(d)),  and  100%  ((e)-(f))  weight  of  spherical  particles. 
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V.  OVERALL  ANISOTROPIC  PROPERTY  CALCULATIONS 


In  this  work,  it  is  assumed  that  the  ellipsoidal  particles  exhibit  transverse  isotropic  be¬ 
havior  {{^e1 :  Ke2 1  Ke3}  =  {400,100,100}  W/111K),  while  the  matrix  and  spherical  particle 
phases  are  isotropic  («:m  =  10  W/mK,  ns  =  50  W/mK).  A  volume  averaged  conductivity 
tensor  for  the  ellipsoidal  phase  is  computed  as 


(23) 


where  R  is  the  rotation  matrix  associated  with  the  orientation  of  the  particle,  and  Xe 
is  the  indicator  function  for  the  ellipsoidal  particles  as  defined  in  Equation  (7).  For  an 
isotropic  distribution  of  transverse  isotropic  ellipsoidal  particles,  the  volume  averaged  tensor 
would  be  isotropic  with  the  mean  thermal  conductivity  Ke  =  200  W/mK.  This  was  verified 
by  considering  uniformly  distributed  rotation  matrices  for  the  given  transverse  isotropic 
conductivity  tensor  [59].  However,  for  an  anisotropic  distribution  of  particles,  the  volume 
averaged  tensor  can  lead  to  a  symmetric  anisotropic  conductivity  tensor. 

As  described  in  Section  III B,  the  integrals  present  in  Equations  (16),  (17)  and  (19) 
are  integrated  numerically  using  the  adaptive  sparse-grid  Smolyak  integration  method  with 
hierarchical  basis.  As  discussed  in  [42] ,  the  numerical  accuracy  depends  on  three  parameters, 
namely  the  quality  of  the  interpolation  defined  as  e,  the  quality  of  the  integration  defined 
as  e,  and  on  the  convergence  of  the  integral,  8%.  The  quality  of  the  interpolation,  e,  is  the 
cutoff  that  determines  which  grid  points  remain  in  the  sparse  grid.  We  use  £  =  lx  10~8 
and  1  x  ICE4  for  the  lower  and  upper  bound  calculations,  respectively.  The  quality  of  the 
integration,  e,  is  related  to  the  volume  removed  in  the  integral  due  to  the  singularity  and  is 


denoted  as 


(24) 


+ 


T f(x  —  x')  (. Srs(x  —  x') 


CrCs)  dfla;/. 


Here,  Ts  and  T f  denote  the  singular  and  formal  terms  of  T,  and  a  discussion  on  the  derivation 
of  these  terms  is  provided  in  Appendix  A.  S(x  —  x')  is  the  Dirac  delta  function,  and  v  is  the 
small  volume  close  to  the  singularity  of  Tf  that  we  remove.  This  volume  v  is  a  sphere  of 
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radius  e  removed  around  the  origin  ( x  =  x').  After  a  convergence  study,  we  selected  e  =  0.1 
mm  for  all  computations.  All  solutions  are  stopped  when  the  integrals  converge  within  a 
certain  tolerance, 

Si  =  \Th  -  lh-i\/\ Lh_x\  <  tol  =  0.001,  (25) 


where  the  subscript  h  is  the  sparse  grid  integration  level  (see  [42]  for  more  details).  In  order 
to  quantify  the  overall  numerical  error  due  to  the  numerical  parameters  mentioned  above, 
we  define 


en  =  x  100[%], 

This  definition  takes  advantage  of  the  property  of  1%,  where 


(26) 


lim  f  Tf(x~x')f(x  —  x')dQxi  =  0  (27) 

V^°J  n.A« 


for  any  statistically  isotropic  function,  f(x  —  x').  In  Equation  (26),  Ksrs(|a._a.'|)  is  computed 
by  numerically  evaluating  the  integral  X  with  isotropic  statistics,  and  ks  is  obtained  in  a 
closed  form  [28].  Considering  the  property  in  Equation  (27),  the  overall  numerical  error  in 
Equation  (26)  measures  the  accuracy  of  the  numerical  integration  of  the  kernel  including  T f 
(see  Equation  (24)).  Note  that  /^s,T.s(i£C_ a,q)  is  not  a  physically  meaningful  quantity  and  is  only 
used  to  determine  the  numerical  error  of  the  calculations.  All  overall  thermal  conductivity 
tensor  calculations  presented  hereafter  are  computed  with  a  numerical  error  en  <  0.16%. 
Note  that  we  are  using  the  material  characteristic  length  scale  to  integrate  over  Vtxi  (see 
discussion  in  sections  IV  A  and  IV  B),  and  a  convergence  study  has  been  performed  to  verify 
the  far  limit  of  the  integration. 

When  computing  the  self-consistent  estimate,  the  NLOPT  library  [60]  (nonlinear  opti¬ 
mization  package)  was  used  for  minimizing  the  objective  functions  (Equation  (19)).  For 
this  work,  the  objective  function  is  considered  converged  when  the  function  value  is  less 
than  1  x  10-6.  Note  that  the  objective  function  is  approaching  0  as  Ko  — )■  77.  It  has  been 
confirmed  that  the  self-consistent  estimates  using  both  the  upper  and  lower  formulations  of 
77  (see  Equation  (19))  yield  identical  results.  This  is  another  verification  test. 

The  resulting  bounds  and  self  consistent  estimate  of  the  overall  anisotropic  thermal  con¬ 
ductivity  tensor  for  the  semi-ordered  pack  of  ellipsoids  are  listed  in  Table  III,  where  the 
eigenvalues  of  the  overall  conductivity  tensor  are  presented.  The  eigenvalues  of  the  con¬ 
ductivity  tensor  (77l,  772,  and  773)  correspond  to  the  components  of  the  conductivity  tensor 
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K\  [W/mK]  «2  [W/mK]  K3  [W/mK]  mean(Ki,K2,K3) 

Isotropic  distribution  of  particles 

Ki  (LB)  - 50.55 — 

Ki  (SC)  - 108.85— 


Ki  (UB)  - 119.47- 

Anisotropic  model 


Ka  (LB) 

59.74 

48.47 

44.11 

50.77 

Ka  (SC) 

141.45 

99.65 

85.27 

108.79 

Ka  (UB) 

154.39 

109.26 

92.23 

118.62 

TABLE  III.  Overall  thermal  conductivity  tensor  predictions  for  semi-ordered  pack. 

in  the  principal  coordinate  frame,  where  the  material  behaves  as  an  orthotropic  one.  Ki 
(subscript  i  stands  for  isotropic)  refers  to  the  assumption  that  the  particles  are  isotropically 
distributed  (perfectly  random),  with  ne  =  200  W/mK,  leading  to  a  macroscopically  isotropic 
conductivity  tensor  with  one  independent  component  ( K\  =  F2  =  k3).  Ki  is  calculated  in  a 
closed  form.  Ka  is  the  anisotropic  conductivity  tensor  (subscript  a  stands  for  anisotropic), 
which  considers  both  the  anisotropic  constituent  and  statistically  anisotropic  configuration. 
Ka  is  calculated  numerically. 

As  can  be  seen  in  Table  III,  the  components  of  the  anisotropic  thermal  conductivity 
tensor  vary  significantly.  For  example,  in  the  upper  bound  of  Ka,  the  difference  between  the 
minimum  and  maximum  eigenvalue  is  154.39  —  92.23  =  62.16  W/mK.  For  the  lower  bound  of 
77a,  we  obtain  59.74  —  44.11  =  15.63  W/mK.  While  the  standard  deviations  in  the  statistical 
functions  were  moderate  (see  Figure  16),  these  deviations  have  a  significant  effect  on  the 
overall  material  behavior  as  they  are  magnified  by  the  anisotropy  in  Ke  (see  Equation  (23)). 
Note  that  the  mean  of  the  eigenvalues  of  Ka  is  nearly  identical  to  results  for  77,;  for  all  three 
models  (UB,  LB  and  SC). 

The  resulting  upper  and  lower  bounds  of  the  thermal  conductivity  tensor  for  the  five 
randomized  packs  and  the  semi-ordered  pack  are  presented  in  Figure  19.  Here,  ^max^min  refers 
to  the  maximum/minimum  eigenvalues  of  the  overall  conductivity  tensor.  These  maximum 
and  minimum  conductivity  components  (solid  lines  with  triangles  for  LB  and  dashed  lines 
with  triangles  for  UB)  are  compared  to  those  assuming  an  isotropic  distribution  of  particles, 
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FIG.  19.  Macroscopic  bounds  for  anisotropic  thermal  conductivity  of  all  compositions  considered. 


Kt  (solid  line  with  circles  for  LB  and  dashed  line  with  circles  for  UB).  Note  that  Kt  (UB)  and 
Kt  (LB)  were  computed  with  k2  =  200  W/mK.  In  general,  the  spread  of  the  bounds  increases 
for  packs  with  more  ellipsoidal  particles,  because  the  difference  between  the  conductivities  of 
an  ellipsoidal  particle  and  the  matrix  is  larger  than  the  difference  between  the  conductivities 
of  a  sphere  and  the  matrix  (see  lines  with  circles  in  Figure  19).  Also  note  that  the  spread 
in  the  minimum  and  maximum  components  of  ~Ra  for  both  bounds  increases  for  packs  with 
more  ellipsoidal  particles  due  to  the  increasing  amount  of  statistical  anisotropy  present  in  the 
packs.  Recall  that  the  statistical  anisotropy  was  quantified  based  on  the  standard  deviations 
in  the  two  point  probability  functions  (see  Figure  18(b),  (d),  and  (f)).  The  spread  in  the 
components  of  the  conductivity  tensor  for  the  lower  bound  is  smaller  than  the  spread  for 
the  upper  bound,  since  the  least  conductive  constituent  (matrix)  is  isotropic  and  the  most 
conductive  constituent  (ellipsoids)  is  anisotropic. 

When  comparing  the  semi-ordered  pack  of  ellipsoids  (filled  markers  on  left  in  Figure  19) 
to  the  randomized  pack  of  ellipsoids,  first  note  that  the  semi-ordered  pack  has  a  larger  Ki 
than  the  randomized  pack,  because  the  semi-ordered  system  has  a  higher  packing  fraction 
(see  subsection  HE 2).  Also  note  that  the  spread  of  the  conductivity  components  is  larger 
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for  the  semi-ordered  pack  due  to  the  larger  statistical  anisotropy  present  in  the  semi-ordered 
pack  (compare  Figures  18(b)  and  16(b)).  We  will  comment  on  the  extent  of  these  differences 
when  we  introduce  a  measure  of  anisotropy  below. 


FIG.  20.  Self-consistent  estimates  for  anisotorpic  thermal  conductivity  of  all  compositions  consid¬ 
ered. 

The  self-consistent  estimates  for  the  mixtures  are  presented  in  Figure  20.  The  same  trends 
described  for  Figure  19  are  observed.  In  order  to  compare  the  self-consistent  estimates  to 
the  LB  and  UB  computations,  we  also  show  the  isotropic  LB  and  UB  behavior,  respectively. 
Note  that  the  self-consistent  estimates  are  closer  to  the  upper  bound  calculations,  since  this 
second-order  model  is  strongly  dependent  on  the  volume  fraction,  and  all  of  the  packs  have 
volume  fractions  greater  than  0.6.  Also  note  that  there  is  no  restriction  on  the  self-consistent 
estimate  being  between  the  lower  and  upper  bounds. 

In  order  to  quantify  the  macroscopic  anisotropy  of  the  thermal  conductivity  tensor,  a 
measure  of  anisotropy,  M^,  is  defined  as 

Ma  =  ma x(K“x  -Ki\,  K™  ~  x  100[%].  (28) 

Computed  measures  of  anisotropy  for  results  in  Figures  19  and  20  are  presented  in  Figure 
21.  As  alluded  to  above,  the  measure  of  anisotropy  increases  as  more  ellipsoidal  particles  are 
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FIG.  21.  Measure  of  anisotropy  for  all  results. 

present  due  to  the  increasing  statistical  anisotropy  in  the  system.  Let  us  consider  the  measure 
of  anisotropy  for  the  upper  bound  (dashed  black  line  in  Figure  21).  Here  =  0.08%  for 
the  100%  mixture  (all  spheres)  and  increases  up  to  =  17.14%  for  the  0%  mixture  (all 
ellipsoids).  For  the  lower  bound,  grows  from  =  0.16%  to  =  5.82%.  The 
measures  of  anisotropy  of  the  upper  and  lower  bounds  for  the  100%  mixture  (all  spheres) 
are  near  the  numerical  error,  \MA(UB )  —  MA(LB)\  =  0.08%  ( en  =  0.16%).  Note  that  the 
anisotropic  upper  and  lower  boimds  are  rigorous  bounds  on  the  overall  anisotropic  thermal 
conductivity  tensor.  Therefore,  these  upper  and  lower  bounds  also  provide  limits  on  the 
measure  of  anisotropy.  For  example,  the  anisotropy  of  the  0%  mixture  (all  ellipsoids)  falls 
between  Myi  =  5.82%  (LB)  and  =  17.14%  (UB). 

The  measures  of  anisotropy  for  the  semi-ordered  pack  (filled  markers  in  Figure  21)  are 
significantly  larger  than  the  estimates  of  the  randomized  packs  and  are  bounded  between 
=  18.19%  (LB)  and  =  29.23%  (UB).  Let  us  consider  the  upper  bound  results,  where 
M^4  =  29.23%  for  the  semi-ordered  pack  (red  filled  square  with  black  outline  in  Figure  21) 
as  compared  to  =  17.14%  for  the  randomized  pack  of  ellipsoids  (dashed  black  line  at 
0%  wt.  spherical  particles  in  Figure  21),  and  let  us  contrast  those  to  std(SVs)  =  0.0072 
for  the  semi-ordered  pack  and  std(5'rs)  =  0.0049  for  the  randomized  pack  (see  section  IV). 
Note  that  the  measure  of  anisotropy  for  the  self-consistent  estimates  of  the  0%  and  25% 
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mixtures  (blue  dotted  dashed  line)  is  higher  than  the  anisotropy  in  the  upper  bound  (black 
dashed  line).  Just  as  a  self-consistent  estimate  is  not  guaranteed  to  be  between  the  upper 
and  lower  bounds,  the  measure  of  anisotropy  associated  with  this  estimate  is  not  expected 
to  be  between  the  upper  and  lower  bounds. 


VI.  CONCLUSIONS 

In  this  work,  we  present  a  systematic  microstructure  characterization  procedure  anchored 
in  micro-CT  data  that  is  used  to  establish  microstructure-statistics-property  relations  of 
polydisperse  particulate  mixtures.  A  novel  image  processing  pipeline  is  developed  that 
accurately  identified  particles  while  maintaining  low  errors.  Improvements  in  the  image 
processing  pipeline  are  achieved  when  compared  to  a  traditional  techniques.  For  all  compo¬ 
sitions  considered,  the  volume  losses  due  to  image  segmentation  are  less  than  4%.  These  low 
errors  indicate  that  scientifically  sound  and  repeatable  results  have  been  achieved.  Next,  we 
developed  a  description  of  the  polydisperse  system  in  continuous  Euclidean  space.  This  ide¬ 
alized  representation  provides  a  substantial  reduction  in  the  dataset  size  and  enables  easier 
data  manipulation  and  understanding. 

After  characterizing  the  microstructure,  n-point  probability  functions  of  real  polydisperse 
mixtures  are  calculated.  We  show  that  second-order  probability  functions  do  not  exhibit 
ellipsoidal  nor  any  other  material  symmetry.  Therefore,  assessment  of  overall  material  con¬ 
stants  in  a  closed  form  is  unattainable. 

The  statistical  description  is  then  used  to  compute  bounds  and  self-consistent  estimates 
of  the  anisotropic  thermal  conductivity  tensor  using  the  Hashin-Shtrikman  variational  prin¬ 
ciple.  This  is  the  first  time  to  the  best  of  our  knowledge,  that  the  anisotropic  second-order 
estimates  are  calculated  without  assumptions  on  the  inclusion’s  shape,  configuration  and/or 
material  anisotropy.  The  overall  properties  show  increasing  anisotropy  in  the  overall  thermal 
conductivity  tensor  for  packs  with  more  transverse  isotropic  ellipsoidal  inclusions.  More¬ 
over,  the  upper  and  lower  bounds  provide  limits  on  the  anisotropy  of  the  mixtures.  Due 
to  the  larger  amounts  of  statistical  anisotropy  for  the  semi-ordered  mixture,  the  measure  of 
anisotropy  for  the  overall  conductivity  tensor  of  this  pack  was  significantly  larger  than  for 
the  randomized  one. 
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Appendix  A:  Discussion  on  r(.x,.x') 


The  infinite  homogeneous  Green’s  function  for  steady-state  heat  conduction  of  an 
anisotropic  material  is 

a(x  ~  x']  =  <d  -  2)1fi(d)fldet  PA  ■  (A1) 

where  det(«)  stands  for  the  determinant,  R  =  ((a;  —  x')  ■  n  ■  (x  —  x'))1^2 ,  d  =  3  is  the 
dimension,  and  f2(d)  is  the  total  solid  angle  in  a  d-dimensional  sphere  [61,  62],  Following  the 
work  of  Torquato  [63],  r(a;  —  x')  is  obtained  by  first  considering  a  temperature  field  defined 
as 

T{x)  =  T0(x)  —  f  V g{x  —  x')  ■  p{x') dfl  (A2) 

Jnx, 

After  integrating  this  equation  by  parts,  using  the  Gauss- divergence  theorem,  and  excluding 
a  small  spherical  region  around  x  =  x',  one  obtains 

Q  —  Qo  +  [  r(aj  —  x')  ■  p(x') dfl  (A3) 

The  resulting  T(aj  —  x')  used  in  Equation  (16)  is  defined  as 

r(a?  —  x')  =  r  f(x  —  x')  —  rs5(a;  —  x'):  (A4) 


where 


’  \x—x'\=Ra 


{g,ihj  +  g./hi)  dfl 


(AS) 


and 


T,(x  -  x') 


d2g 

dxidxj 


(A6) 
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Here,  the  unit  normal  fi  is  the  unit  normal  of  the  spherical  surface  (|£E  —  x'\  =  Rs)  over  which 
the  integration  occurs.  For  an  isotropic  material,  1%  can  be  integrated  analytically,  while 
obtaining  an  analytical  solution  becomes  increasingly  difficult  for  an  anisotropic  one.  In 
this  work,  the  singular  term  is  computed  numerically  using  an  adaptive  multi-dimensional 
integration  scheme  provided  as  an  extension  to  the  GNU  Scientific  Library  (GSL).  The 
singular  term  is  calculated  until  the  relative  error  of  the  integral  is  less  than  1  x  10”5.  Note 
that  the  integral  of  Tf(x  —  x')  over  the  surface  of  a  sphere  is  zero, 

/  Tf(x  -  x')dn  =  0.  (A7) 

J \x—x'\=Rs 


[1]  S.  Kochevets,  T.  J.  J.  Buckmaster,  and  A.  Hegab,  Journal  of  Propellant  and  Power  17,  883 
(2001). 

[2]  K.  Matous,  H.  Inglis,  X.  Gu,  D.  Rypl,  T.  Jackson,  and  P.  Geubelle,  Composites  Science  and 
Technology  67,  1694  (2007). 

[3]  J.  Andrade  and  X.  Tu,  Mechanics  of  Materials  41,  652  (2009). 

[4]  P.  Chaikin  and  T.  Lubensky,  Principles  of  Condensed  Matter  Physics  (Cambridge  University 
Press,  New  York,  2000). 

[5]  J.  Liang  and  K.  Dill,  Biophysical  Journal  81,  751  (2001). 

[6]  P.  Breton,  D.  Ribereau,  C.  Marraud,  and  P.  Lamarque,  International  Journal  of  Energetic 
Materials  and  Chemical  Propulsion  5,  132  (2002). 

[7]  J.  Kepler,  Strena  Sen  de  Niue  Sexangula  (Godefridum  Tampach,  Francofurti  ad  Moenum, 
1611)  translated  as  J.  Kepler,  The  Six-Cornered  Snowflake  (Clarendon,  Oxford,  1966). 

[8]  T.  C.  Hales,  Annals  of  Mathematics  162,  1065  (2005). 

[9]  Y.  Jiao,  F.  Stillinger,  and  S.  Torquato,  Physical  Review  E  79,  041309  (2009). 

[10]  F.  Maggi,  S.  Stafford,  T.  L.  Jackson,  and  J.  Buckmaster,  Physical  Review  E  77,  046107 

(2008). 

[11]  X.  Jia  and  R.  Williams,  Powder  Technology  120,  175  (2001). 

[12]  R.  Caulkin,  X.  Jia,  C.  Xu,  M.  Fairweather,  R.  A.  Williams,  H.  Stitt,  M.  Ni- 

jemeisland,  S.  Aferka,  M.  Crine,  A.  Leonard,  D.  Toye,  and  P.  Marchot, 

Industrial  &  Engineering  Chemistry  Research  48,  202  (2009). 


40 


[13]  S.  Gallier  and  F.  Hiernai'd,  Journal  of  Propulsion  and  Power  24,  147  (2008). 

[14]  B.  Collins,  F.  Maggi,  K.  Matous,  T.  Jackson,  and  J.  Buckmaster,  in  4.6th  AIAA  Aerospace 
Sciences  Meeting  and  Exhibit  (AIAA,  2008)  paper  No:  AIAA  2008-941. 

[15]  T.  Aste,  M.  Saadatfar,  and  T.  Senden,  Physical  Review  E  71,  061302  (2005). 

[16]  H.  Lee,  M.  Brandyberry,  A.  Tudor,  and  K.  Matous,  Physical  Review  E  80,  061301  (2009). 

[17]  C.  Arns,  M.  Knackstedt,  and  W.  Pinczewski,  Geophysics  67,  1396  (2002). 

[18]  J.  Bernal,  Nature  183,  141  (1959). 

[19]  H.  L.  Frisch  and  F.  H.  Stillinger,  Journal  of  Chemical  Physics  38,  2200  (1962). 

[20]  M.  B.  Weissman,  Annual  Review  of  Materials  Science  26,  395  (1996). 

[21]  J.  Schmalzing  and  K.  M.  Grski,  Monthly  Notices  of  the  Royal  Astronomical  Society  297,  355  (1998). 

[22]  C.  Coker  and  S.  Torquato,  Journal  of  Applied  Physics  77,  6087  (1995). 

[23]  C.  Yeong  and  S.  Torquato,  Physical  Review  E  58,  224  (1998). 

[24]  D.  T.  Fullwood,  S.  R.  Niezgoda,  and  S.  R.  Kalidindi,  Acta  Materialia  56,  942  (2008). 

[25]  G.  T.  Seidler,  G.  Martinez,  L.  H.  Seeley,  K.  H.  Kim,  E.  A.  Behne,  S.  Zaranek,  B.  D.  Chapman, 

S.  M.  Heald,  and  D.  L.  Brewe,  Phys.  Rev.  E  62,  8175  (2000). 

[26]  K.  Mecke  and  C.  H.  Arns,  Journal  of  Physics:  Condensed  Matter  17,  S503  (2005). 

[27]  Z.  Hashin  and  S.  Shtrikman,  Journal  of  the  Mechanics  and  Physics  of  Solids  10,  335  (1962). 

[28]  J.  Willis,  Journal  of  the  Mechanics  and  Physics  of  Solids  25,  185  (1977). 

[29]  A.  Bensoussan,  J.  Lions,  and  G.  Papanicolaou,  Asymptotic  Analysis  for  Periodic  Structures. 

(New  York:  North-Holland,  1978). 

[30]  R.  Hill,  Math.  Proc.  Carnb.  Phil.  Soc.  98,  579  (1985). 

[31]  S.  Torquato,  Physical  Review  Letters  79,  681  (1997). 

[32]  R.  Christensen  and  K.  Lo,  Journal  of  the  Mechanics  and  Physics  of  Solids  27,  315  (1979). 

[33]  G.  Dvorak  and  M.  Srinivas,  Journal  of  the  Mechanics  and  Physics  of  Solids  47,  899  (1999). 

[34]  D.  Fullwood,  B.  Adams,  and  S.  Kalidindi,  Journal  of  the  Mechanics  and  Physics  of  Solids 
56,  2287  (2008). 

[35]  L.  Walpole,  Journal  of  the  Mechanics  and  Physics  of  Solids  14,  151  (1966). 

[36]  P.  Areias  and  K.  Matous,  Computer  Methods  in  Applied  Mechanics  and  Engineering  197, 

4702  (2008). 

[37]  C.  Yeong  and  S.  Torquato,  Physical  Review  E  57,  495  (1998). 

[38]  W.  Drugan  and  J.  Willis,  Journal  of  the  Mechanics  and  Physics  of  Solids  44,  492  (1996). 


41 


[39]  S.  Torquato  and  A.  K.  Sen,  Journal  of  Applied  Physics  67,  1145  (1990). 

[40]  P.  Castaneda,  Physical  Review  B  57,  12077  (1998). 

[41]  D.  Talbot  and  J.  Willis,  Journal  of  the  Mechanics  and  Physics  of  Solids  45,  87  (1997). 

[42]  H.  Lee,  A.  Gillman,  and  K.  Matous,  Journal  of  the  Mechanics  and  Physics  of  Solids  56(10), 

1838  (2011). 

[43]  J.  Baker  and  A.  Kudrolli,  Phys.  Rev.  E  82,  061304  (2010). 

[44]  G.  D.  Scott  and  D.  M.  Kilgour,  Journal  of  Physics  D:  Applied  Physics  2,  863  (1969). 

[45]  J.  Russ,  The  Image  Processing  Handbook ,  Image  Processing  Handbook  (Taylor  &  Francis, 

2006). 

[46]  V.  I.  Inc.,  “Amira,”  http://www.amira.com/. 

[47]  B.  microCT,  “Skyscan:  Microtomography,  nanotomography,  non- invasive  3d  xray  mi¬ 
croscopy,”  http : //www . skyscan . be/home . htm. 

[48]  K.  Preston  and  M.  Onoe,  Digital  processing  of  biomedical  images  (Plenum  Press,  1976). 

[49]  G.  Windreich,  N.  Kiryati,  and  G.  Lohmann,  Pattern  Recognition  36,  2531  (2003). 

[50]  J.  Lindblad,  Image  and  Vision  Computing  23,  111  (2005). 

[51]  P.  M.  Chaikin,  A.  Donev,  W.  Man,  F.  H.  Stillinger,  and  S.  Torquato, 

Industrial  <fc  Engineering  Chemistry  Research  45,  6960  (2006). 

[52]  R.  K.  McGEARY,  Journal  of  the  American  Ceramic  Society  44,  513  (1961). 

[53]  A.  Donev,  F.  H.  Stillinger,  P.  M.  Chaikin,  and  S.  Torquato, 

Phys.  Rev.  Lett.  92,  255506  (2004). 

[54]  A.  Donev,  R.  Connelly,  F.  H.  Stillinger,  and  S.  Torquato,  Phys.  Rev.  E  75,  051304  (2007). 

[55]  S.  Torquato,  Random  Heterogeneous  Materials  (Springer,  New  York,  2002). 

[56]  M.  Beran,  Statistical  Continuum  Theories  (Interscience  Publishers,  1968). 

[57]  L.  Onsager,  Phys.  Rev.  37,  405  (1931). 

[58]  G.  Weng,  International  Journal  of  Engineering  Science  30,  83  (1992). 

[59]  J.  Arvo  (Academic  Press  Professional,  Inc.,  San  Diego,  CA,  USA,  1992)  Chap.  Fast  random 
rotation  matrices,  pp.  117-120. 

[60]  S.  G.  Johnson,  “The  nlopt  nonlinear-optimization  package,” 

http : //ab-initio .mit . edu/nlopt. 

[61]  H.-Y.  Kuo  and  T.  Chen,  International  Journal  of  Solids  and  Structures  42,  1111  (2005). 

[62]  Y.  Chang,  C.  Kang,  and  D.  J.  Chen,  International  Journal  of  Heat  and  Mass  Transfer  16,  1905  (1973). 


42 


[63]  S.  Torquato,  Journal  of  the  Mechanics  and  Physics  of  Solids  45,  1421  (1997). 


43 


